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Abstract 

There is compelling evidence that supermassive black holes (SMBHs) exist. Yet the origin 
of these objects, or their seeds, is still unknown. We are performing general relativistic 
simulations of gravitational collapse to black holes in different scenarios to help reveal how 
SMBH seeds might arise in the universe. SMBHs with ~ fO 9 Mq must have formed by 
z > 6, or within 10 9 yrs after the Big Bang, to power quasars, ft may be difficult for gas 
accretion to build up such a SMBH by this time unless the initial seed black hole already has a 
substantial mass. One plausible progenitor of a massive seed black hole is a supermassive star 
(SMS). We have followed the collapse of a SMS to a SMBH by means of 3D hydrodynamic 
simulations in post-Newtonian gravity and axisymmetric simulations in full general relativity. 
The initial SMS of arbitrary mass M in these simulations rotates uniformly at the mass- 
shedding limit and is marginally unstable to radial collapse. The final black hole mass and 
spin are determined to be M^/M « 0.9 and J^/Al^ ~ 0.75. The remaining mass goes into a 
disk of mass M^^/M m 0.1. This disk arises even though the total spin of the progenitor 
star, J/M 2 = 0.97, is safely below the Kerr limit. The collapse generates a mild burst of 
gravitational radiation. Nonaxisymmetric bars or one-armed spirals may arise during the 
quasi-stationary evolution of a SMS, during its collapse, or in the ambient disk about the 
hole, and are potential sources of quasi-periodic waves, detectable by LISA. 



Introduction 

There is substantial evidence that supermassive black holes (SMBHs) of mass ~ 10 6 — 10 10 Mq 
exist and are the engines that power active galactic nuclei (AGNs) and quasars (50J EH HH • 
There is also ample evidence that SMBHs reside at the centers of many, and perhaps most, 
galaxies [55J E5J, including the Milky Way (15112011571 - 

Since quasars have been discovered out to redshift z J> 6 (2H2|, the first SMBHs must have 
formed by zbh ^ 6, or within £bh <, 10 9 yrs after the Big Bang. However, the cosmological 
origin of SMBHs is not known. This issue remains one of the crucial, unresolved aspects of 
structure formation in the early universe. Gravitationally, black holes are strong-field objects 
whose properties are governed by Einstein's theory of relativistic gravitation — general relativity. 
General relativistic simulations of gravitational collapse to black holes therefore may help reveal 
how, when and where SMBHs, or their seeds, form in the universe. Simulating plausible paths by 
which the very first seed black holes may have arisen comprises a timely computational challenge 
(see Fig-QJ. It is an area in which the tools of numerical relativity may be exploited to address a 
fundamental question in cosmology. Performing such simulations may also help identify plausible 
astrophysical scenarios and sites for promising gravitational wave sources involving black holes. 



Figure 1: The formation of a black hole is a strong-held gravitational phenomenon in curved 
spacetime that requires Einstein's equations of general relativity for a description and, in non- 
trivial cases, numerical simulations for a solution. 

We are actively developing new algorithms and new computer codes to solve Einstein's held 
equations of general relativity, coupled to the equations of relativistic hydrodynamics, in three 
spatial dimensions plus time (3+1). As our codes have come online, we have applied them to 
explore different astrophysical scenarios involving strong gravitational helds and the generation 
of gravitational waves. These have included the inspiral and coalescence of binary neutron stars 
and binary black holes, the growth of instabilities in rotating stars, the nonlinear evolution of 
unstable stars, and the collapse and collision of rotating stars and rotating clusters of collisionless 
matter, to name a few. [For recent reviews of simulations of compact binary stars and references, 
see j^ni and We have also performed a wide range of simulations in recent years to study 
alternative scenarios leading to the formation of supermassive black holes, or their seeds. [For a 
review of these calculations and references, see |38| ] Here we will focus on one specihc scenario, 
namely, the collapse of a rotating, supermassive star (SMS) to a supermassive black hole. 

The Onset of Dynamical Instability 

The formation of SMBHs through the growth of black hole seeds by gas accretion is supported 
by the consistency between the total energy density in QSO light and the SMBH mass density in 
local galaxies, adopting a reasonable accretion rest-mass-to-energy conversion efficiency BH] • 
But SMBHs must be present by zbr <; 6 to power quasars. It has been argued [21] that if they 
grew by accretion from smaller seeds, substantial seeds of mass p> 10 5 Mq must already be present 
at z m 9 to have had sufficient time to build up to a typical quasar black hole mass of ~ 10 9 M©. 
A likely progenitor is a very massive object (e.g., an SMS) supported by radiation pressure. 

SMSs (10 3 < M /Mq <^ 10 13 ) may form when contracting or colliding primordial gas builds up 
sufficient radiation pressure to inhibit fragmentation and prevent star formation (see, e.g., |H]). 
SMSs supported by radiation pressure will evolve in a quasi-stationary manner to the point of 



onset of dynamical collapse due to general relativity |51 ll0llT5] Unstable SMSs with M ;> 10 5 Mq 
and metallicity Z <J 0.005 do not disrupt due to thermonuclear explosions during collapse JSj. In 
fact, recent Newtonian simulations suggest that evolved zero- metallicity (Pop III) stars > 300 Mq 
do not disrupt but collapse with negligible mass loss ^7]. This finding could be important since 
the first generation of stars may form in the range 10 2 — 10 3 Mq [51 Q] . A combination of turbulent 
viscosity and magnetic fields likely will keep a spinning SMS in uniform rotation f|H H71l5T?ll3^| : 
but see [23 and the final section below for an alternative). As they cool and contract, uniformly 
rotating SMSs reach the maximally rotating mass-shedding limit and subsequently evolve in a 
quasi-stationary manner along a mass-shedding sequence until reaching the instability point. At 
mass-shedding, the matter at the equator moves in a circular geodesic with a velocity equal to 
the local Kepler velocity 2 . 

It is straightforward to understand the radial instability induced by general relativity in a 
SMS by using an energy variational principle |5()1 I42| . Let E = E{p c ) be the total energy of a 
momentarily static, spherical fluid configuration characterized by central mass density p c . The 
condition that E(p c ) be an extremum for variations that keep the total rest mass and specific 
entropy distribution fixed is equivalent to the condition of hydrostatic equilibrium and establishes 
the relation between the equilibrium mass and central density: 



dE 
dp c 



= M cq = M cq (p c ) (equilibrium). (1) 



The condition that the second variation of E(p c ) be zero is the criterion for the onset of dynamical 
instability. This criterion shows that the turning point on a curve of equilibrium mass vs. central 
density marks the transition from stability to instability: 

d 2 E dM ca 

—^r = <=> — = (onset of instability). 2 

dpc 2 op c 

Consider the simplest case of a spherical Newtonian SMS supported solely by radiation pres- 
sure and endowed with zero rotation. This is an n = 3, (T = 1 + 1/n = 4/3) polytrope, with 
pressure 

P = Prad = ~aT 4 = Xpt, (3) 

where K = K(s la d) is a constant determined by the value of the (constant) specific entropy 
Srad = ^aT 3 /n in the star. Here T is the temperature, n is the baryon number density, and a is 
the radiation constant. Consider a sequence of configurations with the same specific entropy but 
different values of central density. The total energy of each configuration is 

E{ Pc ) = E/ rad + W, (4) 

where J7 ra d is the total internal radiation energy and W is the gravitational potential energy. 
Applying the equilibrium condition Q to this functional yields M cq = M cq (s ra d), i.e. the equi- 
librium mass depends only on the specific entropy and is independent of central density (see 
Fig- Eh)- Applying the stability condition then shows that all equilibrium models along this 
sequence are marginally stable to collapse. 

Now let us account for the effects of general relativity. If we include the small (de-stabilizing) 
Post-Newtonian (PN) correction to the gravitational field, we must also include a comparable 
(stabilizing) correction to the equation of state arising from thermal gas pressure: 

P = frad + P e ,s = ^aT 4 + 2nkT 7 (5) 
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Figure 2: A sketch of mass versus central density along an equilibrium sequence of SMSs of 
fixed entropy. Panel (o) shows nonrotating, spherical Newtonian models supported by pure 
radiation pressure; (b) shows nonrotating, spherical PN models supported by radiation pressure 
plus thermal gas pressure; (c) shows rotating PPN models spinning at the mass-shedding limit. 



where we have taken the gas to be pure ionized hydrogen. Note that P ga8 /-Prad = 8/(s ra d/fc) -c 1. 
The energy functional of a star now becomes 

E{p c ) = U, ad + W + AU gils + AW PN , (6) 

where AU gas is the internal energy perturbation due to thermal gas energy and AWpn is the PN 
perturbation to the gravitational potential energy. Applying the equilibrium condition Q now 
yields M oq w M^ wt times a slowly varying function of p c (see Fig.[2J)). The turning point on the 
equilibrium curve marks the onset of radial instability; the marginally stable critical configuration 
is characterized by 

p c , crit = 2 x 10- 3 M 6 " 7/2 gm cm" 3 , 

T c , crit = (3 x lO^Mg- 1 K, (7) 
(R/M) crit = 1.6 X 10 3 M 6 1/2 , 

where M 6 denotes the mass in units of 10 6 M Q . (Here and throughout we adopt gravitational 
units and set G = 1 = c.) 

Finally, let us consider a uniformly rotating SMS spinning at the mass-shedding limit [2]. A 
centrally condensed object like an n — 3 polytrope can only support a small amount of rotation 
before matter flys off at the equator. At the mass-shedding limit, the ratio of rotational kinetic 
to gravitational potential energy is only T/|W| = 0.899 x 10~ 2 <C 1. Most of the mass resides in 
a nearly spherical interior core, while the low-mass (Roche) envelope bulges out in the equator: 
i? oq /i?poi c = 3/2. When we include the contribution of rotational kinetic energy to the energy 
functional, we must now also include the effects of relativistic gravity to Post-Post-Newtonian 
(i.e. PPN) order, since both T and AWpn scale with p c to the same power. The energy functional 
becomes 

E{ Pc ) = U Iad + W + A[/ gas + AP^pn + AT^ppn + T (8) 

Applying the equilibrium condition lfl]l. holding M , angular momentum J and s fixed, now yields 
M cq ps M^ ewt times a slowly varying function of p c (see Fig. Et)- If we restrict our attention 
to rapidly rotating stars with M > 10 5 Mq (the typical size of the seed black holes adopted in 
some recent galaxy merger simulations; see, e..g., [22]) the influence of thermal gas pressure is 
unimportant in determining the critical point of instability. The turning point on the equilibrium 



curve then shifts to higher density and compaction than the critical values for nonrotating stars, 
reflecting the stabilizing role of rotation: 




Pc,crit 



0.9 x 10 _1 Mg" 2 gm cm 
(9 x 10 7 )M 6 " 1/2 K, 



-3 



(9) 



The actual values quoted above for the critical configuration were determined by a careful numer- 
ical integration of the general relativistic equilibrium equations for rotating stars !2-; they are in 
close agreement with those determined analytically by the variational treatment. The numbers 
found for the nondimensional critical compaction and angular momentum are quite interesting. 
First, they are universal ratios that are independent of the mass of the SMS. This means that 
a single relativistic simulation will suffice to track the collapse of a marginally unstable, maxi- 
mally rotating SMS of arbitrary mass. Second, the large value of the critical radius shows that 
a marginally unstable configuration is nearly Newtonian at the onset of collapse. Third, the fact 
that the angular momentum parameter of the critical configuration J/M 2 is below unity suggests 
that, in principle, the entire mass and angular momentum of the configuration could collapse to 
a rotating black hole in vacuum without violating the Kerr limit for black hole spin (but see the 
next section below!). 

The Outcome of Collapse 

There are several plausible outcomes that one might envision a priori for the dynamical collapse 
of a uniformly rotating SMS once it reaches the marginally unstable critical point identified 
above. It could collapse to a clumpy, nearly axisymmetric disk, similar to the one arising in the 
Newtonian SPH simulations for the isothermal (r = 1) implosion of an initially homogeneous, 
uniformly rotating, low-entropy cloud |25| . Alternatively, the disk might develop a large-scale, 
nonaxisymmetric bar. After all, the onset of a dynamically unstable bar mode in a spinning 
equilibrium star occurs when the ratio T/|W| ~ 0.27 (see, e.g., JI] and |2J for Newtonian 
treatments and fjlH] and for simulations in general relativity). Since T/|W| is 0.899 x 10~ 2 at 
the onset of collapse and scales roughly as i?" 1 during collapse, assuming conservation of mass 
and angular momentum, this ratio climbs above the dynamical bar instability threshold when 
the SMS collapses to R/M « 20, well before the horizon is reached. The growth of a bar might 
begin at this point. Indeed, a weak bar forms in simulations of rotating supernova core collapse 
|281 17] . but in supernova the equation of state stiffens (r > 4/3) at the end of the collapse, 
triggering a bounce and thereby allowing more time for the bar to develop. A rapidly rotating 
unstable SMS might not form a disk at all, but instead collapse entirely to a Kerr black hole; not 
surprisingly, a nonrotating spherical SMS has been shown to collapse to a Schwarzschild black 
hole . Alternatively, the unstable rotating SMS might collapse to a rotating black hole and 
an ambient disk. 

Recently we have performed two simulations that together resolve the fate of a marginally 
unstable, maximally rotating SMS of arbitrary mass M. In (SEj we followed the collapse in full 
3D, but assumed PN theory. We tracked the implosion up to the point at which the central 
spacetime metric begins to deviate appreciably from flat space at the stellar center. We found 
that the massive core collapses homologously during the Newtonian epoch of collapse, and that 
axisymmctry is preserved up to the termination of the integrations. This calculation motivated us 
45 to follow the collapse in full general relativity by assuming axisymmctry from the beginning 
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Figure 3: Snapshots of density and velocity profiles during the implosion of a marginally unstable 
SMS of arbitrary mass M rotating uniformly at break-up speed at t = 0. The contours are drawn 
for p/p m ax = lCP ' 43 0=0 — 15), where p max denotes the maximum density at each time. The 
fourth figure is the magnification of the third one in the central region: the thick solid curve at 
r « 0.3M denotes the location of the apparent horizon of the emerging SMBH. (From Shibata & 
Shapiro gS]) 



to maximize spatial resolution (see Fig- ■ We found that the final object is a Kerr-like black 
hole surrounded by a disk of orbiting gaseous debris. The final black hole mass and spin were 
determined to be M^/M w 0.9 and J/JM^ w 0.75. The remaining mass goes into the disk of mass 
Ma^/M « 0.1. A disk forms even though the total spin of the progenitor star is safely below the 
Kerr limit. This outcome results from the fact that the dense inner core collapses homologously 
to form a central black hole, while the diffuse outer envelope avoids capture because of its high 
angular momentum. Specifically, in the outermost shells, the angular momentum per unit mass 
j, which is strictly conserved on cylinders, exceeds jisco, the specific angular momentum at 
the innermost stable circular orbit about the final hole. This fact suggests how the final black 
hole and disk parameters can be calculated analytically from the initial SMS density and angular 
momentum distribution |40j . The result applies to the collapse of any marginally unstable n = 3 
polytrope at mass-shedding. Maximally rotating stars which are characterized by stiffer equations 
of state and smaller n (higher T) do not form disks, typically, since they are more compact and 
less centrally condensed at the onset of collapse 02] ■ 

The above calculations show that a SMBH formed from the collapse of a maximally rotating 
SMS is always born with a "ready-made" accretion disk. This disk might provide a convenient 
source of fuel to power the central engine. The calculations also show that the SMBH will be 
born rapidly rotating. This fact is intriguing in light of suggestions that observed SMBHs rotate 
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Figure 4: Quasi-stationary evolution snapshots of a SMS that cools and contracts in the absence 
of any viscosity or magnetic fields from a nearly spherical initial state that is in very slow, uniform 
rotation. The maximum density is normalized to unity; the highest density contour level is 0.9 and 
subsequent levels range from 10 _1 to 10~ 10 and are separated by a decade. The final configuration 
has reached the mass-shedding limit and is toroidal and differentially rotating, with T/|W| « 0.27 
and n polc /fl cq fa 2.5 x 10 3 . (From New & Shapiro [77] 



rapidly (e.g., gElE]). 

The implosion will result in a burst of gravitational waves. Current relativistic collapse cal- 
culations for this problem break down soon after the apparent horizon of the hole approaches 
its asymptotic size. Moreover, the need to resolve the centrally condensed inner core precludes 
extending the numerical grid very far out into the wave zone. These two restrictions prevent a 
reliable determination of the complete burst waveform at present. Crude estimates can be gener- 
ated for the frequency and amplitude based on the quadrupole formula, and the results suggest 
that they are well within the range of detectability for LISA for a reasonable spectrum of masses. 
More interesting, perhaps, is the possible generation of quasi-periodic gravitational waves from 
nonaxisymmetric instabilities that might arise in the ambient disk; the characteristics of these 
waves are likely to be comparable those discussed in the next section and in Eqn. 1|10|) below. 

An Alternative Scenario 

The previous analysis assumes that a combination of turbulent viscosity and magnetic fields 
combine to keep a spinning SMS rotating uniformly prior to the onset of collapse. While this is 
the most likely possibility El IH3 EHI i New and Shapiro |27] have discussed an alternative. In 
the absence of viscosity or magnetic fields, an evolving SMS will conserve angular momentum on 
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Figure 5: Snapshots of density contours in the equatorial plan of differentially rotating Newtonian 
toroids of different polytropic indicies, n. All of the stars are constructed from the same differential 
rotation law with r2 po i e /f2 eq ss 26 and T/|W| ~ 0.14 at t = 0. The snapshots are taken after the 
stars have evolved for many central rotation periods, P c . Stars with higher n exhibit a one-armed 
spiral m=l instability. The panels show results for (n,t/P c ) — (a)(2, 37); (6)(2.5, 24); (c)(3, 17); 
and (d)(3.33, 19). The contour lines denote densities p/p c = 10"( 16 ^ d (t = 1, • • • , 15). (From 
Saijo, Baumgarte & Shapiro 

cylinders as it cools and contracts in a quasi-stationary manner. So even if it is rotating uniformly 
at formation, it will evolve to a state of differential rotation. Our evolution calculations indicate 
that an axisymmetric configuration that is slowly rotating uniformly and nearly spherical at birth 
(T/|1K| <C 1) will evolve to a differentially rotating toroid (see Fig.QJ). Prior to encountering any 
relativistic radial instability, the star reaches the mass-shedding limit, at which point T/|W| ps 
0.27 and differential rotation is extreme, f£ po ie/^eq ~ 2.5 x 10 3 . 

Toroidal configurations which have high values of T/|W| and substantial differential rotation 
and central mass concentration are subject to nonaxisymmetric dynamical instabilities, especially 
the one-armed spiral m = 1 instability Softer equations of state with higher polytropic 

indicies are very susceptible to this instability (see Fig. [SJ, so this mode is likely to be triggered 
in a rotating SMS should it evolve with low viscosity and magnetic field. Such a mode typically 
induces a secondary m=2 bar mode of smaller amplitude, and the bar mode can excite quasi- 
periodic gravitational waves (see Fig. El . Typical wave amplitudes and wave frequencies for a 
source at a distance r are given by 

h ~ 10- 16 (T/|^|)(M/i?)A/ 6 r G i c ^ 10- 20 Af 6 r G i c , (10) 
/ Gw ~ lQ- l {T/\W\) l/2 {M/Rf /2 Ms l Hz~l(r A Ms l Hz. 

where we have set T/| W| ~ 0.27 and R/M ~ 10 3 in evaluating the above quantities. These values 
also apply to waves generated by bar modes that may arise in a similar fashion in the ambient 
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Figure 6: Gravitational radiation waveforms as seen by a distant observer on the rotation axis for 
the configurations shown in Fig. |S] The radius r is the distance to the source and the ratio M/ R 
is the (arbitrary) compaction of the initial configuration. (From Saijo, Baumgarte & Shapiro |34|) 

disk surrounding the SMS black hole that forms from the collapse of a relativistically unstable, 
uniformly rotating SMS (see discussion above). Both the amplitudes and frequencies put these 
waves well within the detectable range for LISA for a realistic distribution of masses. Future 
simulations that we are planning will explore this possibility in greater detail, and help resolve 
which scenario best describes the final fate of a rotating SMS. 
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